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Abstract 



Using operator algebra, we extend the series for the activity density 
in a one-dimensional stochastic sandpile with fixed particle density 
p, the first terms of which were obtained via perturbation theory [R. 
Dickman and R. Vidigal, J. Phys. A 35, 7269 (2002)]. The expansion 
is in powers of the time; the coefficients are polynomials in p. We 
devise an algorithm for evaluating expectations of operator products 
and extend the series to 0(t 16 ). Constructing Pade approximants to a 
suitably transformed series, we obtain predictions for the activity that 
compare well against simulations, in the supercritical regime. 
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I. INTRODUCTION 



Sandpiles with a strictly conserved particle density (so-called fixed- energy sand- 
piles or FES [1]), exhibit absorbing-state phase transitions [2-4], and have recently 
attracted much interest. Until now, most quantitative results for FES have been 
based on simulations [5-8], an important exception being the solution by Priezzhev 
et al. [9] of a directed, fixed-energy version of the Maslov-Zhang model [10], via the 
Bethe ansatz. It is therefore of great interest to develop theoretical approaches for 
FES. 

Series analysis has proved to be one of the most accurate and reliable approaches 
to critical phenomena, in both equilibrium and nonequilibrium contexts [11-14]. 
Series expansion typically functions best in low- dimensional systems (because longer 
series can be derived), that is, for just those systems in which the renormalization 
group and expansion about an upper critical dimension d c are less reliable. In the 
case of sandpiles, systematic epsilon expansions are as yet unavailable, and the 
value of d c in fact remains controversial [15-17]. Simulation results suggest novel 
critical behavior in the one-dimensional FES, although conflicting critical exponent 
values have been reported [7,8,18,19], which may reflect finite-size effects. Series 
expansions, on the other hand, implicity treat the infinite-size limit, and so provide 
important information, complementary to that afforded by simulations. In light of 
these observations, we believe it highly desirable to apply series methods to sandpile 
models. 

This paper is one of a series analyzing a stochastic sandpile using operator meth- 
ods. In an earlier work [20] a path-integral representation was developed and an ex- 
pansion derived for the order parameter (activity density) in powers of time. While 
the path-integral formalism reveals interesting features of the model, and may be 
applied in any number of dimensions, the complexity of the diagrammatic expansion 
limits the number of terms that can be obtained. (In Ref. [20] terms up to 0(t 5 ) 
are reported.) In this paper we employ a different approach, which permits us to 
extend the series for the one- dimensional case considerably. After casting the master 
equation for the sandpile in terms of an operator formalism, we analyze the direct 
expansion of its (formal) solution, leading to an algorithm for generating the series 
coefficients. 

We consider Manna's stochastic sandpile in its fixed-energy (particle-conserving) 
version [7,20-22]. The configuration is specified by the occupation number n at each 
site; sites with n > 2 are said to be active, and have a positive rate of toppling. When 
a site topples, it loses exactly two particles ("grains of sand"), which move randomly 
and independently to nearest-neighbor (NN) sites. (Any configuration devoid of 
active sites is absorbing, i.e., no futher evolution of the system is possible once such 
a configuration is reached.) In this work, as in Ref. [20], we adopt a toppling rate 
of n(n— 1) at a site having n particles, which leads us to define the order parameter 
as p = {n{n — 1)). While this choice of rate represents a slight departure from the 
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usual definition (in which all active sites have the same toppling rate), it leads to 
a much simpler evolution operator, and should yield the same scaling properties 
[20]. Preliminary simulation results [23] indicate that in one dimension the model 
exhibits a continuous phase transition at p c = 0.9493. 

In the following section we define the model and review the operator formalism 
introduced in Ref. [20]. This is followed in Sec. Ill by an analysis leading to 
an expansion in terms of so-called reduced commutators. Implementation of the 
expansion in a computational algorithm is described in Sec. IV. Then in Sec. V 
we report numerical results of the series analysis. A summary and discussion is 
provided in Sec. VI. 



II. MODEL 



As discussed in Ref. [20], the master equation for this model may be written in 
the form 



m 

dt 



L|M>), 



(1) 



where 



I*) = £p(KM)IM>, 

{n;} 



is the probability distribution, and the evolution operator takes the form, 



-(7Ti_i + n i+1 ) 2 - 7lf 



(2) 



Here a« and rii are, respectively, annihilation and creation operators associated with 
site i, defined via 

aAm) = nAni-l) 



and 



7TiK) = K + l) 



The formal solution of the master equation is \^(t)) = e tL \^/(0)); that for the activity 
density is: 



Pit) = ( K^mo)). 



(3) 
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Here we have introduced the notation: 

<I=£<WI (4) 
M 

for the projection onto all possible states; thus normalization reads: ( |^) = 1. We 
consider a uniform Poisson-product initial distribution. Letting p n = e~ p p n /n\, and 
using \P)i = J2 ni P ni \ n )i to denote a Poisson distribution at site i, we have, 

i*(o)> =ni p >*- ( 5 ) 



We shall expand equation (3) for the activity density in powers of t. 



III. OPERATOR ALGEBRA 

To begin we note some basic properties of operators a,j, iTj and Lf 



a r ;\P) 3 =p n \P) 3 (6) 

( K- = ( I (7) 

{\Li = (8) 

[<H,^j] = S id (9) 

[Li, Lj] = [<n, Lj] = for \i - j\ > 1 . (10) 



The second relation expresses the fact that the creation operator conserves the 
normalization of any state, while the third shows that Lj conserves probability, as 
it must. 

The coefficient of t n /n\ in the expansion of the activity is: 

J2i \a%L Sl L S2 ...L Sn \P) , (11) 
s 

where the sum is over all sequences S of sites s = 0,«i,,,s n with \s±\ < 1, and 
s j+1 E {s jtTnin - 1, s j)Tnax + 1}, for j > 1, where s j)min = min{s , ...Sj}, and s jtmax 
is the maximum of this set. The restriction on sequences follows from equations (8) 
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and (10); if the condition were violated, it would be possible to move one of the Lj 
to the left of all other operators, yielding a result of zero. 

Our strategy for evaluating p s is to commute each Lj to the left of eft. The 
first step replaces ciqL s1 by its commutator, due to equation (8). If we write this 
commutator in normal order, that is, with all creation operators 7ij to the left of all 
annihilation operators, then the 7r's may be replaced by 1, by equation (7). Thus, 

(\alL 3 = (\[alL 3 ] R , (12) 



where the subscript R denotes a reduced commutator, that is, the commutator in 
normal order, with all 7r's replaced by unity. Evidently [oq, Lj] R involves only anni- 
hilation operators. The two nontrivial expressions of this kind are: 

[a 2 , L ] R = -2a 2 - 4aj* , (13) 



and 



[al,Lx} R = -a{ + 2a a\ 



(14) 



In the computational algorithm, discussed in some detail below, it is not neces- 
sary to generate the tree structure of sequences explicitly, since each monomial is 
processed separately and both translation and reflection symmetries may be used in 
the calculations of the contributions ps- Evaluating the expectation of each term in 
ps is trivial, because 

(|a-.<"|P)=p M , (15) 



where M = ^ . rrij is the number of annihilation operators, irrespective of which 
sites are involved. 

It remains to find a general expression for the reduced commutator [F^\ Lk] R . 
bmce is a linear combination of products of annihilation operators, and recalling 
that <2j and Lk commute if \i—k\ > 1, we see that the problem reduces to evaluating 

C(p, q, r) = L ] R . (16) 



(Commutators involving Lj with j j^O are obtained using translation invariance.) It 
is straightforward to evaluate C(p, q, r) using the following identities. First we note 
that 

[a p j ,7i j }=pa p r 1 , (17) 
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as is readily shown by induction. Using this it is simple to show 



[dj,Hj]R = p(p-l)dj 2 + 2pa P j 1 



18) 



and 



[ag,nM] R =p(p-lX + 2pa* 



(19) 



Finally, we may use equation (17) to show that for i^j, 



(20) 



Applying these relations one readily finds: 
C(p,q,r)=al +2 



\p(p-l)arM + \r(r-iy_ ia r* 



■ ^ p— 1 r—1 i p— 1 r i P r— 1 

+ -W a -\ a i + pd-i OL\ +ra_ 1 a\ 
- qa^al [2a^ +1 + (q-l)al\ . 



(21) 



Using this result, we can evaluate the reduced commutators in a computer algorithm. 



IV. COMPUTATIONAL ALGORITHM 



Let us discuss some details of the computer algorithm used to generate the series 
for the activity. We employ a recursive procedure to generate the contributions 
of order n + 1 on the basis of those of order n. From equation (3) and with a 
Poisson-product initial distribution defined in equation (5) we notice that 

% = ( \alLe tL \P) = (\[al,L] R e tL \P) = -p+( \(Aa a\ - 4ajj) e tL \P) . (22) 



The last equality above may be understood using the reduced commutators (13) 
and (14). Using reflection symmetry, we have ( | [oq, L] r e tL \P) = ( |{[oq, L ] r + 
2[al, Li] R }e tL \P) and further simplification is provided by translation symmetry. 
The coefficient of t n /n\ in the expansion of the activity may be identified with C n , 
where 



C n — 



d n p 



dt r 



(23) 



t=Q 



Using the procedure described above we obtain a recursion relation 
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Cn+1 — — C n + (|F n+ i|P), 



(24) 



where F n+1 = [F n , L] and iq = Aaoa 2 — 4oq. To complete the algorithm, we have that 
Co(0) = p(0) = ( |ag|) = p 2 , where relation (15) is used. An immediate consequence 
of these recursion relations is that the coefficient of p 2 in the term of order n in 
t is given by (—1)™, since each monomial in the functions F has at least three 
annihilation operators [24]. This was already shown in Ref. [20]. 

The calculations were done in two steps. Initially, the functions F were calcu- 
lated up to order 12. Each monomial was represented by three integer variables: 
an eight-byte integer for the numerator of the coefficient, a four- byte integer for the 
denominator (which is always a power of 2), and another eight-byte integer to store 
the number of factors of each annihilation operator. Since all calculations are done 
in integer arithmetic, there are no roundoff errors. Using translation invariance, 
each monomial was put in a form such that the the annihilation operator of lowest 
spatial index is a . The power mi associated with each annihilation operator a Si (as 
in equation (15)) is stored in four bits of the eight-byte integer variable mentioned 
above. As each new monomial is generated, a search is performed for any existing 
term with the same set of powers; storing all powers in a single integer facilitates 
the search. As a consequence of equation (10), when the reduced commutator of a 
monomial with L is calculated, nonzero contributions may arise only from the com- 
mutator of the monomial with L_i, L , Li, . . . , L i+ i, where % is the largest index in 
the monomial. As is clear from equation (21), each of these commutators can give 
rise to up to seven new monomials. Thus, it is apparent that the number of mono- 
mials grows very rapidly as the order is increased; the function F 12 involves 519 115 
monomials. To go beyond order 12 it is necessary to handle integers larger than 
can be represented using eight bytes and to process monomials with more than 16 
exponents, which can no longer be stored in a single 8-bite integer variable. In fact, 
at order 12 most of the processing time is used in the search procedure. Therefore, 
our results from order 13 to 16 were obtained processing the monomials in F 12 one- 
by-one, generating all contributions from it at orders 13-16. In these calculations 
the numerators were represented by two eight-byte integers. The limiting order (16) 
is determined by the large number of new monomials generated; a single monomial 
in F15 may generate on the order of 40 monomials in F^. The results presented here 
required about 170 hours of cpu time on an Athlon K7 1800 MHz computer. 

It is convenient to write the expansion in the form: 

««) s # = E^Ew- (25) 

n m=0 



The series coefficients b 1hm are listed in Table I. In Ref. [20] it was shown that 
m < n — 1, with 



(2«-l)!! 



*•■"- = 2 WTW: (26) 
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The coefficients reported in Table I satisfy this relation at each order, and agree 
with those derived (for n < 5) using the path-integral formalism [20]. 



V. ANALYSIS OF SERIES 



The coefficients b n ^ m jn\ in the time series, equation (25), grow rapidly with n; 
the rate of growth appears to be faster than exponential. This is evident from an 
analysis of 

Km = ^ ) (27) 



To see if the h rhm follow a systematic trend we analyze these quantities for a given 
q = (m — l)/(n — 2). (For a fixed value of n, the h n ^ m appear to trace out a smooth 
curve, so that h n (q) for intermediate values of q can be estimated via interpolation.) 
As shown in figure 1, h n (q) appears to grow faster than exponentially with n, away 
from the limits q = and q — 1. (Observe that for m — 0, h m , n — > — oo as 
n — > oo since b n $ = 1, and similarly for m — n — 1, since equation (26) implies that 
b n ,n-i grows more slowly than n\ .) A reasonable description of the dominant growth 
in the series coefficientes is h n (q) ~ n a (q), with the exponent a (see the inset of 
figure 1) taking its maximum value of about 1.2 for q ~ 0.4. This of course implies 
faster-than-exponential growth for the coefficients b n ^ m /n\. 

Next we examine the behavior of the coefficients in the time series for specific 
values of the particle density p. Let 

?(f) = 5>H) B (28) 

n 



where 

n-l 



= E b ^P m ■ (29) 



m=0 



For p = 1 (slightly above the critical value of 0.9493), c n is simply the sum of all 
coefficients at order n, divided by n\. The coefficients c n again grow faster than 
exponentially, with lnc n ~ n 1 ' 15 for p — 1 and ~ n 1 - 10 for p — 2. (Given the limited 
number of coefficients, we cannot make very precise estimates of the exponent. The 
key point is that the growth appears to be faster than exponential.) These results 
imply that equation (28) is a divergent series with zero radius of convergence. 

We turn now to an analysis of the series for p(t). As is well known, it is often 
possible to obtain useful results from divergent series by means of a resummation 
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technique. In the present case we construct Pade approximants to the time series or 
to the series obtained via a transformation of variable [13,26-28]. We have examined 
many transformations, for example 

l-e~ bt , x 

y = , (so) 



x = — > — , (31) 
1 + bt K ' 



z=l- t, — T~r , ( 32 ) 
(1 + btp ' v ; 



w = l t-, rr , (33) 

1 + ln(l + bt)' K ' 



and 



v = 1 -exp{6[l - (1 + ty<]} . (34) 



Each transformation maps the interval t > to a finite interval, and can be 
expanded as a power series in t about t — 0, with the lowest-order term oc t. Each is 
readily inverted permitting one to express the time t in powers of the new variable. 
The slow convergence associated with the power-law or logarithmic forms in the 
last four expressions is motivated by the numerical finding of slow relaxation in the 
sandpile model, even far from the critical point [25]. We analyze the transformed 
series for p (or for hip) using Pade approximants. The degree of success depends 
greatly on the range of p under consideration. (Each transformation includes a free 
parameter b, which can be adjusted to optimize the regularity of the result, or to 
obtain consistency between different approximants. Except where noted, the results 
do not exhibit much sensitivity to the choice of this parameter.) 

For small values of p, the best results are obtained via Pade approximants to the 
t-series without any transformation of variable. Figure 2 compares series predictions 
for p = 0.5 (obtained using the [6,7], [7,8] and [7,9] approximants to the series for 
p(t)) against the result of a Monte Carlo simulation for a system of 500 sites (for 
p = 0.5 finite-size effects are negligible at this system size). The [7,8] approximant 
is reliable for t < 10. (Various other approximants, such as [8,8] and [7,7], are 
ill-behaved and provide reasonable predictions only for quite short times, typically 
t < 2.) We have not been able to improve the series prediction for longer times, 
either by a change of variable or through analysis of lnp(t) or its time derivative. 
Although some improvement could be expected with longer series, it appears unlikely 
that the asymptotic decay of p(t) in the subcritical regime will be accessible through 
analysis of an expansion in powers of time. 
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For larger values of p the transformation defined in equation (32), using 7 = 1/2, 
is the most useful of those studied. In figure 3 we compare the [8,8] approximant 
(obtained using b = 0.57) with simulation data for p — 1. The situation is markedly 
better than for p — 1/2: the series prediction accompanies the simulation result 
up to around t = 1000. It must be noted, however, that the good agreement seen 
here depends on the choice of the transformation parameter b. For other values the 
agreement with simulation is not as good. (A more suitable criterion for choosing 
b would be by seeking agreement among various approximants [13]. In the present 
case this is not possible because the off-diagonal approximants to the z series are ill- 
behaved, while the [7,7] approximant behaves very similarly to the [8,8] used here.) 
Despite the good agreement up to times on the order of 1000, the present series 
seems incapable of capturing the asymptotic long-time relaxation of p{t), which is 
non-monotonic, as shown in figure 3. 

Remarkably, the reliability of the series improves dramatically at larger values of 
the particle density p. Series and simulation results for ~p(t) at p = 2 are compared 
in figure 4; the maximum relative error is about 0.1%. (The series prediction is gen- 
erated as for p — 1, but using b = 1.5 in this case.) The good agreement, moreover, 
persists at long times, motivating a study of p^ = lim t ^ 00 ~p{t) , corresponding to the 
transformed series with z — lva. equation (32). (We again use the [8,8] approximant 
to the ^-series.) 

The series prediction for p^ is compared with simulation in figure 5, using pa- 
rameters b = 0.57 and 6 = 5. Excellent agreement is found for p > 2, the relative 
error being < 0.2%. The smaller b value yields better results for p ~ 1, whereas 
slightly better results are obtained for larger p, using the larger b value. For p > 2 
we may claim quantitative accuracy for the series prediction. Nearer the critical 
point, the agreement appears reasonable (at least on the scale of figure 5), but it is 
clear that the 16-term series cannot be used to study critical properties. For exam- 
ple, the prediction using b = 0.57 yields a critical value of about 0.906, that is, the 
extrapolated activity density goes to zero at this p value. The critical value found 
in simulations is 0.9493. 

In summary, the present series seems quite reliable in the supercritical regime, 
both at short and at asymptotically long times, whereas its utility in the critical 
and subcritical regimes is restricted to rather short times. Just above the critical 
point, rather good predictions are possible for short and intermediate times, but 
this depends on a judicious choice of the transformation parameter b. 



VI. DISCUSSION 

We develop an algebraic method leading to a time series for the activity density 
of the stochastic sandpile model introduced in [20]. Determination of the series 
coefficients depends on evaluation of certain commutators, an algebraic task readily 



10 



codified in a computational algorithm. We extend the series for the one-dimensional 
case to sixteen terms. 

Analysis of the series yields disappointing results for the subcritical and critical 
regimes, but very good predictions in the supercritical region, as judged by com- 
parison with Monte Carlo simulation. At first glance this is surprising, since in the 
subcritical regime the stationary state is inactive and might be regarded as trivial. 
Relaxation to this inactive state (and to the active state at or near the critical point 
p c ) is however nontrivial, characterized by stretched exponential, power-law, or other 
slowly-converging forms [25] . It appears to be very difficult to capture such behavior 
in the kind of temporal series developed here, which employs a Poissonian initial 
distribution. The reason is that for smaller values of the particle density (p < 2, 
say), the one-site stationary occupation distribution P(n) is far from Poissonian. As 
p increases, the second factorial moment (n(n — 1)) = p approaches p 2 , as expected 
for a Poisson distribution. Figure 7 shows that the stationary one-site distribution 
observed in simulations approaches the corresponding Poisson distribution with the 
same density p. (Even for p = 8 there are significant differences between the distri- 
butions; but analysis of the third and fourth factorial moments suggests convergence 
to a Poisson distribution as p — > oo.) 

An important open question is whether simply increasing the number of terms 
would permit one to analyze the small-p regime. The present results suggest that 
even with 20 or 30 terms this region would remain inaccessible. It appears to be 
more promising to approach the critical region from above, since for larger particle 
densities we find good agreement with numerical results. In this context it is in- 
teresting that the quality of predictions near the critical point improves greatly on 
going from 12 to 16 terms. This suggests that further extension of the series, to 20 
or more terms, would yield quantitative results for critical properties, through study 
of relaxational properties at p c , or or of stationary properties (as in figure 5) as p c is 
approached from above. It would also be of great interest to develop an expansion 
for a stationary property such as directly in powers of the particle density p, but 
this appears to be much more difficult than deriving an expansion in powers of the 
time. We leave such investigations, using modified or extended series, as subjects 
for future work. 
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TABLES 



n 


m 


bn,m 








1 


1 





1 


2 





1 




1 


8 


3 





1 




1 


66 




2 


80 


4 





1 




1 


442 




2 


2 076 




3 


896 


5 





1 




1 


2 842 




2 


35 396 




3 


52 240 




4 


10 752 


6 





1 




1 


18 118 




2 


516 880 




3 


1 737 952 




4 


1 187 968 




5 


135 168 


7 





1 




1 


922 265 
8 




2 


7 040 282 




3 


45 847 512 




4 


67 368 480 




5 


25 614 368 




6 


1 757 184 


8 





1 




1 


5 865 473 
8 




2 


370 752 137 
4 




3 


1 078 168 434 




4 


2 871 388 040 




5 


2 283 464 832 




6 


536 472 640 




7 


23 429 120 


9 





1 




1 


74 596 747 






16 




2 


4 797 745 191 
4 




3 


23 841 662 132 




4 


105 679 404 154 



13 





5 


147 137 780 760 




6 


71 353 965 088 




7 

8 


11 072 770 560 
318 636 032 


10 





1 




1 

2 


474 336 627 




16 

123 077 063 429 




3 


8 

1 018 938 641 745 
2 




4 


3 584 915 570 625 




5 


7 999 349 570 432 




6 


6 656 488 808 368 




7 

8 


2 121 777 710 528 
227 436 059 136 




9 


4 402 970 624 


11 





1 




1 

2 
3 


12 064 410 263 




64 

3 142 928 518 289 
16 

85 419 503 179 415 

8 

116 020 128 091 449 




4 




5 


394 806 480 115 048 




6 


514 548 057 479 072 




7 

8 


278 154 455 793 952 
61 313 513 593 600 




9 


4 683 285 856 256 




10 


61 641 588 736 


12 





1 




1 

2 
3 


76 711 895 439 




64 

80 070 040 225 479 
32 

1 770 456 755 814 995 




4 


8 

14 610 068 149 248 089 
4 




5 


18 396 700 126 638 476 




6 


35 650 110 284 461 928 




7 

8 


29 745 976 515 005 712 
11 054 665 928 232 448 




9 


1 747 506 609 502 464 




10 


97 252 577 107 968 




11 


872 465 563 648 
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3 
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2 037 418 656 354 491 
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16 

905 058 014 398 112 835 
8 

1 655 391 460 208 555 433 
2 




6 


2 309 179 626 832 648 726 




7 


2 816 714 002 502 804 952 



14 





8 


1 601 275 099 838 022 656 




9 


426 223 203 786 122 496 




10 


49 655 626 778 919 936 




11 


2 046 635 410 882 560 




12 


12 463 793 766 400 
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1 




1 

2 
3 


2 613 736 799 297 
64 

11 934 019 637 184 639 
32 

711 799 150 376 749 517 




4 


8 

53 725 847 102 644 113 051 
16 

142 451 880 202 934 178 839 
4 




5 




6 


140 971 191 603 315 396 510 




7 

8 


243 993 717 303 561 711 492 
201 186 302 850 192 322 944 




9 


81 746 928 038 823 569 408 




10 


16 113 035 142 846 829 824 




11 


1 415 730 263 534 155 776 




12 


43 811 063 460 921 344 




13 


179 478 630 236 160 


15 





1 




1 

2 
3 
4 


531 822 407 449 409 
2048 

606 748 047 325 325 193 
128 

116 603 968 592 784 196 927 
64 

819 559 563 865 455 675 379 




5 


8 

12 371 190 775 573 187 034 899 

8 




6 


8 499 377 166 784 889 887 638 




7 

8 


20 286 352 375 993 324 998 496 
23 298 464 567 721 533 566 328 




9 


13 579 307 980 015 469 945 184 




10 


4 068 154 005 082 098 401 408 




11 


607 110 051 667 479 652 352 




12 


40 856 912 571 394 580 480 




13 


957 525 442 027 462 656 




14 


2 602 440 138 424 320 
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1 




1 

2 
3 


422 699 161 810 361 
256 

61 687 281 835 997 869 817 
1024 

1 192 545 870 991 479 699 681 




4 


32 

12 465 559 773 385 531 628 949 




5 


4 

133 184 916 649 036 025 384 179 
2 




6 


502 866 857 598 243 404 511 546 




7 

8 


1 627 359 034 988 855 536 002 199 

2 537 319 446 210 036 202 445 148 




9 


2 040 132 355 769 944 077 918 400 




10 


872 811 268 389 569 306 302 976 
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12 


22 853 161 358 227 259 040 768 




13 


1 195 547 596 367 062 589 440 




14 


21 401 594 847 260 721 152 




15 


37 965 009 078 190 080 



Table I. Series coefficients in the expansion of the activity. 
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FIGURE CAPTIONS 



Figure 1. Function h n (q) as defined in text, for q = 0.2 (■); q = 0.4 (•); q = 0.6 

(o); q — 0.8 (□). Observe that h grows faster than linearly for q = 0.2, 0.4. Inset: 
growth exponent a(q) defined via h n (q) ~ n a( - q \ 

Figure 2. Normalized activity p = p(t)/p 2 versus time for p — 1/2 from simulation 
and various Pade approximants to the times series as indicated. 

Figure 3. Normalized activity for p = 1. Symbols: simulation result; curve: series 
prediction as described in text. 

Figure 4. As in figure 3 for but for p = 2. 

Figure 5. Main graph: limiting activity versus particle density p. Points: simu- 
lation; solid curve: series prediction using transformation (32) with b = 0.57, [8,8] 
Pade approximant; dashed line: same approximant and transformation but using 
6 = 5. Inset: difference A = ~p^ series - for b = 0.57 (■) and b = 5 (□). 

Figure 6. Single-site occupancy distributions P{n) obtained in simulation (■) com- 
pared with the corresponding Poisson distribution (□). Upper panel: p = 1.2; 
middle: p = 3; lower: p = 8. 
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